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Executive  summary 


/ - \ 

Background 

Metamaterials  and  smart  active  structures  offer  a  huge  potential  to  develop  multifunctional  systems 
addressing  the  technological  needs  of  army  forces.  Like  evolutionary-optimised  biological  systems  found  in 
Nature  (e.g.  adaptive  skin  texture  of  cephalopods),  active  soft  structures  producing  large  deformations  offer 
attractive  ways  to  incorporate  the  ability  to  reconfigure  their  surface  texture.  Surfaces  that  modify  their 
(potentially  multiscale)  texture  in  response  to  either  direct  control  or  automatic  feedback  from  external  cues 
(solar  light  intensity,  temperature,  incoming  waves,  etc)  are  particularly  attractive  as  they  could  be  used  to 
change  radar,  acoustic,  thermal  or  optical  signatures  to  fulfil  mission-specific  objectives. 


Objectives 

A  natural  question  would  be:  how  does  one  design  a  surface  that  can  dynamically  change  shape  and  have  its 
texture  smoothly  and  reversibly  reconfigured  into  specified  complex  three-dimensional  patterns  at  specific 
spatial  scales?  How  to  control  the  deformation  of  active  elements  made  of  metamaterials  to  obtain  these 
textured  patterns? 

The  goal  of  this  project  was  to  provide  a  proof-of-concept  answer  to  this  question  and,  by  doing  so,  effectively 
bridging  the  gap  between  the  integration  of  metamaterials  and  the  design  of  reconfigurable  active  structures. 


Approach 

Here,  it  was  proposed  to  develop  a  robust  computational  modelling  platform  based  on  non-linear  finite 
element  and  optimisation  procedures  to  design  and  optimise  patterned  surfaces  by  using  a  clever  combination 
of  materials,  structural  layers,  boundary  and  loading  conditions  (active  controls).  The  idea  was  to  control 
three-dimensional  surfaces  by  exploiting  their  local  buckling/wrinkling/folding  characteristics  triggered  by 
smooth  deformation  of  active  reinforcement  elements. 


Results 

The  proof-of-concept  approach  was  successfully  demonstrated.  It  is  possible  to  deform  the  surface  of  a  multi¬ 
layer  soft  structure  with  a  high  level  of  control  to  match  a  given  target  three-dimensional  surface  topography 
by  inducing  controlled  localised  buckling.  The  multiphysics  computational  platform  that  was  developed  is 
general  and  can  accommodate  various  types  of  active  materials  acting  as  controllable  actuators. 

s _ * 


Keywords :  Camouflage,  Biomimetics,  Active,  Adaptive,  Morphing,  Surface,  Structure,  Multi-Scale  Structural 
Mechanics,  Metamaterials,  Theoretical  and  Computational  Methods,  Finite  Element,  Isogeometric,  NURBS, 
Optimisation 
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1.  Introduction 


1.1  Background 

At  the  micro-  and  nanoscale,  the  periodicity  of  patterned  surfaces  present  very  interesting  properties  such  as 
unique  microfluidic  [1],  optical,  electronic  or  acoustic  properties  because  of  wave  interference  and/or 
transformation  [2].  In  Nature,  a  number  of  animals  and  plants  have  evolved  the  ability  to  exploit  patterned 
surfaces  at  various  scales  (e.g.  Lotus  flower,  shark  and  its  skin  cuticles)  to  their  advantages  [3],  Dynamic  or 
adaptive  patterned  surfaces  are  very  desirable  as  they  can  be  adapted  to  specific  contexts,  operating 
conditions  and  environments.  Cephalopods  such  as  octopus  (Figure  1),  squid  and  cuttlefish  possess  some  of 
the  most  striking  abilities  in  the  animal  kingdom:  camouflage  via  colour  change  and  adaptive  change  of  the 
three-dimensional  texture  of  their  skin  [4-6].  In  a  military  context,  these  characteristics  are  very  attractive  and 
can  open  up  a  vast  array  of  applications. 


Figure  1.  Octopus  demonstrating  its 
camouflage  abilities  in  its  natural 
environment  (Sea  of  Cortez,  Mexico,  G. 
Limbert). 


1.2  Military  relevance 

Metamaterials  and  smart  active  structures  are  particularly  relevant  to  the  technology  requirements  of  army 
forces.  The  ever  growing  trend  is  to  develop  autonomous  (also  self-healing  and  fault-tolerant)  systems  that  can 
provide  multifunctionality  whilst  minimising  weight/size,  particularly  in  extreme  environments  (underwater, 
air  and  space).  The  solution  to  these  challenges  is  to  develop  single  systems  that  can  be  reconfigured  for 
specific  missions  /  applications.  Like  evolutionary-optimised  biological  systems  found  in  Nature  (e.g. 
cephalopods),  active  soft  structures  producing  large  deformations  offer  attractive  ways  to  incorporate  the 
ability  to  reconfigure  their  surface  texture.  Surfaces  that  modify  their  texture  in  response  to  either  direct 
(human/computer-controlled)  orders  or  automatic  feedback  from  external  cues  (solar  light  intensity, 
temperature,  incoming  waves,  etc)  would  provide  the  army  forces  with  a  distinct  technological  advantage  in  a 
wide  variety  of  scenarios  and  applications.  This  would  make  such  multifunctional  active  structures  particularly 
attractive  as  they  could  be  used  to  change  radar,  acoustic,  thermal  or  optical  signatures  to  fulfil  specific 
mission  requirements. 

A  natural  question  would  be  how  does  one  design  a  surface  that  can  dynamically  change  shape  and 
have  its  texture  smoothly  and  reversibly  reconfigured  into  specified  complex  three-dimensional  patterns 
at  specific  spatial  scales?  How  to  control  the  deformation  of  active  elements  made  of  metamaterials  to 
obtain  these  textured  patterns ? 

The  goal  of  this  project  was  to  provide  a  proof-of-concept  answer  to  this  question. 

Here,  it  was  proposed  to  develop  a  robust  computational  modelling  platform  based  on  non-linear  finite 
element  and  optimisation  procedures  to  design  and  optimise  patterned  surfaces  by  using  a  clever 
combination  of  materials,  structural  layers,  boundary  and  loading  conditions  (active  controls).  The  idea  was  to 
control  three-dimensional  surfaces  by  exploiting  their  local  buckling/wrinkling/folding  characteristics  triggered 
by  smooth  deformation  of  active  reinforcement  elements.  Dynamic  reconfiguration  of  the  texture  could  be 
done  by  using  electroactive  [7]  or  magnetoactive  materials  [8], 
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The  proposed  proof-of-concept  innovation  could  be  used  to  assist  engineers  in  DOD-related  industries  and 
laboratories  to  accelerate  the  understanding,  design  and  development  of  smart  active  surfaces  by  exploiting, 
in  a  robust  and  systematic  way,  the  physics  of  metamaterials  [e.g.  electroactive  polymers  (EAPs)]  to  propose 
innovative  solutions  for  military  forces.  This  would  effectively  close  the  gap  between  the  integration  of 
metamaterials  and  the  design  of  reconfigurable  active  structures.  The  main  unique  selling  point  of  the  concept 
proposed  in  this  project  is  to  provide  and  integrated  computational  modelling  environment  so  that  novel  cost- 
effective  solutions  could  be  explored  in  a  minimum  amount  of  time  (computational  optimisation  studies). 

1.3  Proposed  approach 

The  proposed  computational  platform  is  based  on  non-linear  finite  element  and  numerical  optimisation 
techniques. 

Because  of  the  novelty  and  highly  non-linear  nature  of  the  problems  to  be  simulated  (see  Remark  1),  new 
and/or  improved  finite  element  techniques  and  numerical  algorithms  had  to  be  developed  for  this  project 
(they  are  detailed  in  section  2).  To  efficiently  capture  wrinkling  instabilities,  it  was  proposed  to  develop  a 
robust  isogeometric  (NURBS-based)  finite  element  framework  [9]  (work  package  1,  section  2.2)  and  an 
improved  arc-length  non-linear  finite  element  solver  (work  package  2,  section  2.3).  A  special  optimisation 
procedure  had  to  be  developed  to  establish  local  active  deformation  patterns  of  the  top  surface  of  a  multi¬ 
layer  structure  required  to  produce  a  final  target  shape  (work  package  3,  section  2.4).  For  subsequent 
exploitations  of  the  computational  framework  post-processing  and  visualisations  algorithms  were  also 
developed  (work  package  2,  section  2.5). 

Finally,  the  computational  platform  was  exploited  and  proof-of-concept  results  are  provided  in  section  3. 

Remark  1 

Mechanics  of  wrinkle  formation:  the  theoretical  study  of  wrinkling  has  recently  received  much  attention 
in  the  physics  and  soft  matter  communities  [10-13]  as  they  are  relevant  for  a  wide  range  of  industrial 
applications  and  fundamental  problems.  Wrinkles  form  in  response  to  competing  physical  constraints 
such  as  inextensibility,  mismatch  of  stiffness  in  multi-layer  structures  and  specific  multi-axial  loading 
conditions  together  with  a  principle  of  energy  minimisation  [12].  Wrinkling  of  initially  flat  (linearly) 
elastic  thin  sheets  and  their  non-linear  kinematics  are  governed  by  the  fourth-order  non-linear  Foppl-von 
Karman  equations  [14]  which  cannot  be  solved  analytically  in  most  cases,  particularly  if  one  considers 
complex  geometry,  boundary  and  loading  conditions.  The  bending  energy  term,  critical  to  wrinkling,  is 
proportional  to  the  second  derivative  of  the  local  curvature  of  the  shell  mid-surface  and,  in  a  finite 
element  context,  this  means  that  C1- continuity  across  finite  element  boundaries  is  required  [15, 16], 


Remark  2 

In  the  original  project  proposal,  it  was  proposed  to  conduct  a  high-throughput  computer  experiment  to 
assess  the  characteristics  of  wrinkles  created  by  using  a  large  number  of  combinations  of  material  and 
structural  parameters  for  a  bi-layer  structure  as  well  as  various  uniform  and  non-uniform  boundary  and 
loading  conditions.  From  this  design  of  computer  experiment,  a  surrogate  model  (also  called  meta¬ 
model)  would  have  been  built  and  used  to  calculate  the  sensitivities  of  the  physical  system  to  the  design 
variables  as  well  as  cross  correlations  between  these  variables  [17],  It  was  hypothesised  that  the 
exploitation  of  the  surrogate  model  would  provide  fundamental  insights  into  the  relations  existing 
between  wrinkle  patterns,  wavelength,  amplitude,  mechanical  and  underlying  structural  properties.  By 
analysing  the  direct  sensitivities  and  output  response  of  the  surrogate  model  predictive  laws  would  have 
been  established.  The  surrogate  model  would  have  been  subsequently  used  to  conduct  (computationally 
cheap)  multiobjective  numerical  optimisation  to  seek  specific  wrinkle  characteristics  (and  therefore 
specific  three-dimensional  texture)  [17, 18). 

This  idea  was  later  replaced  by  that  of  conducing  direct  optimisation  studies  to  determine  parameters 
necessary  to  obtained  particular  target  surfaces  (see  sections  2.4  and  2.5). 
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2.  Software  development 


For  this  project,  extensive  software  development  was  conducted  and  consisted  of  the  following  main  work 
packages: 

■  WP1:  Development  of  isogeometric  finite  element  formulations  and  their  numerical  implementation 

■  WP2:  Development  of  inverse  non-linear  finite  element  optimisation  procedures 

■  WP3:  Development  and  implementation  of  a  robust  arc-length  solver  for  wrinkling  analysis 

■  WP4:  Development  and  implementation  of  advanced  algorithms  for  tensor  field  visualisation 

The  computational  platform  was  exploited  to  address  the  objectives  of  the  research  project,  namely 
demonstrating  the  feasibility  of  controlling  the  3D  texture  of  multi-layered  soft  structures  by  inducing  localised 
buckling  behaviour  (see  section  3) 

2.1  General  software  development  platform  for  mathematical  formulations  and  numerical 
implementations 

In  the  initial  phase  of  the  project,  the  implementation  of  an  isogeometric  thin  shell  element  and  associated 
finite  element  procedures  were  carried  out  using  the  high-level  language  and  numerical  environment  of 
MATLAB®  (MathWorks,  Inc.,  Natick,  MA,  USA)  (see  Appendix  A  -  Journal  paper  submitted).  The  choice  of  this 
particular  platform  was  mainly  made  because  of  the  extensive  experience  of  the  post-doctoral  researcher  in 
using  MATLAB®. 

In  the  second  phase  of  the  project,  a  comparable  environment  based  on  Mathematica®  (Wolfram  Research, 
Inc.,  Champaign,  IL,  USA)  was  thoroughly  used.  This  environment  allowed  us  to  significantly  accelerate  the 
development  and  this  choice  will  prove  very  judicious  indeed  for  any  future  iterations  of  the  project.  A 
symbolic-numeric  Mathematica®  plugin  suite  AceGen/AceFEM  (http://www.fgg.uni-lj.si/Symech/)  was  used. 
The  Mathematica®  package  AceGen  is  used  for  the  automatic  derivation  of  formulae  needed  in  numerical 
procedures.  Symbolic  derivation  of  the  characteristic  quantities  (e.g.  gradients,  tangent  operators,  sensitivity 
vectors,  etc)  generally  leads  to  exponential  growth  of  derived  expressions,  both  in  time  and  space.  A  very 
robust  and  efficient  approach  implemented  in  AceGen  avoids  this  problem  by  combining  several  techniques: 
symbolic  and  algebraic  capabilities  of  Mathematica®,  automatic  differentiation  technique,  automatic  code 
generation,  simultaneous  optimisation  of  expressions  and  theorem  proving  by  stochastic  evaluation  of  the 
mathematical  expressions  [19].  The  multi-language  capabilities  of  AceGen  can  be  used  for  a  rapid  prototyping 
of  numerical  procedures  in  script  languages  of  general  problem  solving  environments  like  Mathematica  or 
MATLAB®  as  well  as  to  generate  highly  optimised  and  efficient  compiled  language  codes  in  Fortran  or  C. 
Through  a  unique  user  interface  the  derived  formulae  can  be  explored  and  analysed. 

The  AceGen  package  also  provides  a  collection  of  prearranged  modules  for  the  automatic  creation  of  the 
interface  between  the  automatically  generated  code  and  the  numerical  environment  where  the  code  would  be 
executed.  The  AceGen  package  directly  supports  several  numerical  environments  such  as:  MathLink 
connection  to  Mathematica®,  AceFEM  which  is  a  research  finite  element  environment  based  on  Mathematica® 
or  research  and  commercial  finite  element  packages  based  on  the  Fortran  language.  The  multi-language  and 
multi-environment  capabilities  of  the  AceGen  package  enable  the  generation  of  numerical  codes  for  various 
numerical  environments  from  the  same  symbolic  description.  In  combination  with  the  finite  element 
environment  AceFEM  the  AceGen  package  represents  an  ideal  tool  for  rapid  development  of  new  numerical 
formulations  and  algorithms.  The  visualisation  platform  that  was  used  was  ParaView  (www.paraview.org). 


\ 

HH ParaView 

Visualisation 

J 


Figure  2.  Software  environment  used. 
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2.2  WP1:  Development  of  isogeometric  finite  element  formulations  and  implementation 

This  work  package  consisted  in  implementing  a  generic  modular  framework  for  2D  and  3D  isogeometric  finite 
element  formulations  [9]  into  the  environment  of  AceFEM  (see  section  2.1).  As  mentioned  in  the  Introduction 
section,  mechanical  instabilities  involving  wrinkling  and  folding  feature  high  curvature  gradients  that  can  only 
be  fundamentally  captured  by  using  C1  formulations  [16]  and  these  imply  continuity  of  the  displacement  field 
across  boundaries  of  finite  elements.  Most  "traditional"  Lagrange  polynomial-based  finite  element 
formulations  do  not  satisfy  these  requirements  unlike  elements  based  on  B-Spline  and/or  non-uniform  rational 
B-Splines  (NURBS)  [9,  20],  The  idea  of  using  NURBS  [21]  as  basis  functions  for  analysis  was  introduced  by 
Hughes  et  al.  [9,  20],  and  was  named  isogeometric  analysis.  In  isogeometric  analysis,  the  functions  for  the 
geometry  description  are  also  used  as  basis  functions  for  the  analysis.  Thus,  the  analysis  works  on  a 
geometrically  exact  model.  This  offers  the  possibility  to  close  the  existing  gap  between  design  and  analysis  by 
merging  design  geometry  and  analysis  model.  It  was  demonstrated  that  not  only  were  NURBS  applicable  to 
engineering  analysis,  but  that  they  were  better  suited  for  many  applications,  and  were  able  to  deliver  accuracy 
superior  to  standard  finite  elements  (see,  e.g.,  [22-30]).  More  importantly,  NURBS  are  smooth,  higher  order 
functions  which  have  become  standard  in  CAD  (computer  aided  design)  programs.  They  allow  great  geometric 
flexibility  and  high  order  continuities  at  the  same  time.  NURBS  are  therefore  ideally  suited  as  basis  functions 
for  Kirchhoff-Love  shell  [15,  31-33], 

To  be  compatible  with  traditional  data  structures  used  in  finite  element  software  applications,  the  architecture 
of  the  isogeometric  numerical  implementation  was  based  on  the  use  of  a  Bezier  extraction  operator  [34]  which 
basically  maps  NURBS  to  Bezier  polynomial  bases  and  therefore  offers  a  "traditional"  finite  element  data 
structure  for  elements. 


The  developed  isogeometric  framework  works  for  NURBS-based  geometries  defined  by  tensor  products  up  to 
cubic  polynomial  order.  The  framework,  which  is  a  collection  of  specially  written  Mathematica®  functions  and 
C  codes,  can  handle  2D  and  3D  multiphysics  (i.e.  multi-fields)  element  formulations  and  can  accommodate 
material  constitutive  laws  of  arbitrary  complexity.  Natural,  essential  and  initial  boundary  conditions  are 
seamlessly  integrated  without  the  need  to  use  penalty  functions  to  enforce  them.  Edge  and  surface  loads  are 
also  accounted  for.  By  making  use  of  the  powerful  differentiation  capabilities  of  AceGen,  all  algorithms  and 
element  formulations  are  consistently  linearised  so  that  non-linear  solution  procedures  converge  quadratically 
[35], 

Dedicated  direct  sensitivity  subroutines  were  also  developed.  In  a  finite  element  context,  the  object  of  a 
sensitivity  analysis  is  to  calculate  the  dependence  of  a  functional  response  5R(q)  =  {^>1(q),^,.,(q),...,^,(q)} 
made  of  s  scalars  on  n  parameters  q  =  {qvq2,...,q  }  which  are  generally  arbitrary  analysis  model  inputs 
(constitutive  parameters,  geometric  characteristics  of  the  finite  element  mesh,  boundary  and  loading 
conditions)  but  can  also  be  arbitrary  intermediate  or  final  results  of  the  analysis.  This  is  achieved  by  calculating 
the  local  partial  derivatives  of  the  functional  response  with  respect  to  the  parameters  leading  to  a  sensitivity 
matrix  &  : 


6  =  [1] 

<9q 

The  physical  interpretation  of  a  calculated  sensitivity  is  that  it  provides  the  quantitative  change  of  the  response 
variable  as  result  of  a  unit  change  in  the  parameter  considered.  Sensitivities  can  therefore  be  exploited  in 
numerical  optimisation  studies  and  accelerate  convergence  of  the  procedure  (this  approach  is  followed  in 
WP2,  section  2.4).  In  the  context  of  this  project  numerical  subroutines  were  developed  to  calculate  three  types 
of  sensitivity: 

•  Primal  sensitivity  analysis:  sensitivity  of  the  primal  field  variables  (3  displacement  degrees  of 
freedom  for  each  node  of  the  mesh  (u,  v,  w) )  with  respect  to  the  constitutive  parameters. 

•  Shape  sensitivity  analysis:  sensitivity  of  the  primal  field  variables  (3  displacement  degrees  of 
freedom  for  each  node  of  the  mesh  (u,v,w))  with  respect  to  the  geometrical  characteristics  of 
the  mesh  (and  the  location  of  the  boundary  conditions). 

•  Secondary  sensitivity  analysis:  sensitivity  of  selected  response  variables  (e.g.  Cauchy  stress 
tensor  [a  ,a  ,cr  .a  ,  <J  and  &  ]  and  von  Mises  stress)  with  respect  to  the  constitutive 

L  XX  '  yy  '  ZZ  '  xy  '  XZ  yz  J  1 

parameters. 
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In  this  work  package,  post-processing  Mathematica®  functions  were  also  developed  to  project  NURBS  element 
results  to  linear  2D  and  3D  elements  so  that  results  of  isogeometric  analyses  could  be  visualised  in  ParaView 
(see  section  2.5).  A  dedicated  mean-square  fitting  of  Gauss-points  fields  with  user-controlled  mesh  refinement 
was  devised. 


Figure  3.  Deformed  shape  of  a  thin 
structure  after  a  highly  non-linear 
isogeometric  analysis  post-buckling 
procedure.  Colour  plot  corresponds 
to  the  first  principal  Green- 
Lagrange  strains  projected  from  the 
NURBS  elements  of  the 
isogeometric  analysis  to  a  mesh 
consisting  of  tri-linear  hexahedron 
elements. 


Deliverables 

■  NURBS-based  Isogeometric  finite  element  framework  for  primal  and  direct  sensitivity  analyses,  and 
post-processing  for  static  and  dynamic  non-linear  solving  procedures  based  on  implicit  scheme  [35], 

■  2D  and  3D  continuum  elements  implemented  (linear,  quadratic  and  cubic  polynomial  order)  with 
possibility  to  couple  kinematically  2D  and  3D  elements. 

■  Structural  element:  Mindlin-Reissner  general  shell  which  degenerates  to  a  Kirchhoff-Love  thin  to 
ultra-thin  shell  in  the  limit  of  vanishing  shell  thickness  without  the  typical  associated  pathological 
locking. 

■  Periodic  boundary  conditions  implemented. 


2.3  WP3:  Development  and  implementation  of  a  robust  arc-length  solver  for  wrinkling  analysis 

Post-buckling  behaviour  of  structures  (e.g.  wrinkling)  is  a  typical  type  of  non-linear  behaviour  where  a  standard 
Newton-Raphson  procedure  can  fail.  Arc-length  solution  methods  offer  a  good  alternative  [36,  37],  An  arc- 
length  iterative  solution  procedure  solves  consistently  a  linearised  augmented  system  which  is  composed  of 
equilibrium  equations: 

R(Ap,AA)  =  0  m 


and  constraint  equation: 


ff(Ap,  AA)  =  0 


[3] 


for  a  set  of  variables  p  and  A  with  unknown  values  which  respectively  correspond  to  the  vector  containing 
the  global  nodal  degrees  of  freedom  and  load  multiplier.  The  constraint  equation  is  defined  as: 


g( Ap,AA)  =  ApTAp  +  AA2^|R^,/R.re/  +  AA2'P^p'Ic/prc/  -  As  =  0  [4] 

and  considers  an  increment  of  arc-length  of  the  equilibrium  path  As  to  determine  the  incremental  variables 
Ap  and  AA . 

rfj  and  \l/2  are  scaling  parameters  while  R rej  and  prej  are  respectively  the  reference  load  vector  (natural 
boundary  conditions)  and  the  reference  vector  of  prescribed  displacements  (essential  boundary  conditions). 


Generally,  arc-length  methods  are  based  on  a  two-step  solution  strategy.  The  first  step  is  a  predictor  phase 
which  executes  the  first  iteration  and  defines  the  direction  of  the  path  following  procedure  (positive  or 
negative)  in  the  current  increment.  For  this  reason,  a  judicious  criterion  to  determine  the  correct  direction  is 
required.  The  second  step  is  a  corrector  phase  and  executes  the  subsequent  iterations  in  the  current 
increment. 

The  implementation  of  the  arc-length  solver  into  AceFEM  was  based  on  the  work  of  de  Souza  Neto  and  Feng 
[38]  for  an  efficient  criterion  to  determine  the  path  direction  and  followed  a  similar  approach  to  that  of 
Schweizerhof  and  Wriggers  [39].  Consistent  linearisation  of  the  solution  procedure  ensures  quadratic  rate  of 
convergence  of  the  Newton-type  algorithm. 
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2.4  WP2:  Development  of  inverse  non-linear  finite  element  optimisation  procedures 

2.4.1.  Rationale  and  description 

The  optimisation  procedure  was  developed  to  establish  what  local  activation  control  patterns  of  embedded 
active  elements  (embedded  as  an  array  of  reinforcement  bars  within  a  soft  polymeric  structure  made  of  one  or 
several  layers)  are  necessary  to  deform  a  surface  to  match  a  specified  3D  shape. 

Based  on  the  choice  of  material  properties  and  geometric  characteristics  of  a  multilayer  soft  surface  containing 
particular  arrangements  of  active  contractile/extensible  elements  (that  could  be  made  of  electroactive 
polymers),  the  user  can  choose  a  3D  target  shape  that  is  specified  by  providing  an  analytical  formula  but  which 
could  equally  be  obtained  by  reverse  engineering  real  surfaces  using  laser  or  tomographic  techniques.  The 
computational  platform  is  then  used  to  run  an  inverse  finite  element  optimisation  study  that  calculates  the 
localised  active  deformations  required  in  each  reinforcement  bar  element  to  produce  the  targeted  shape. 


2.4.2.  Active  elements 

For  this  proof-of-concept  project,  the  real  physical  phenomena  driving  the  active  linear  contraction/extension 
of  the  reinforcement  elements  were  not  important  to  model.  Unidimensional  deformations  of  the 
reinforcement  bar  elements  were  prescribed  by  specifying  a  single  contraction/extension  scalar  parameter  (j> 
so  that  the  stretch  of  the  element  was  AA  =  An  +  <f>  with  A0  being  the  initial  stretch  values  (equal  to  1  in  a 
stress-free  configuration).  This  is  based  on  a  simple  thermal  analogy.  Elements  were  implemented  with 
quadratic  shape  functions  and  used  Lobatto  integration  rules  [35,  36], 

Kinematic  coupling  procedures  were  developed  to  embed  bar  elements  (by  projection)  into  shell  elements 
composing  the  top  surface  of  the  multi-layer  structure.  A  special  projection  technique  was  implemented  so 
that  any  type  of  reinforcement  bar  element  geometry  could  be  used. 


2.4.3.  Numerical  optimisation 

In  this  project,  only  linear  constraints  were  considered  for  the  optimisation  problem: 

■  Minimum  and  maximum  elongations  of  reinforcement  bar  elements 

■  Targeted  displacement  field  required  to  morph  the  original  multi-layer  structure  into  the  targeted 
shape 

■  Maximum  deformation  of  the  mesh  of  each  solid  layer. 

The  objective  function  and  all  constraints  are  linear  functions  of  the  optimisation  variables  (</r),  i  being  the 
identifier  of  each  reinforcement  bar  element).  The  problem  is  therefore  a  linear  programming  (LP)  problem. 
This  type  of  optimisation  problem  can  be  efficiently  solved  using  a  simplex  method  [40],  Solution  procedures 
for  non-linear  optimisation  problems  can  be  viewed  a  sequence  of  linear  problems  and  can  be  amenable  to  a 
LP  form.  The  optimum  solution  of  an  LP  problem  corresponds  to  one  of  the  basic  feasible  solutions  and  thus 
can  be  found  by  examining  all  basic  solutions.  However,  depending  on  the  particular  problem,  the  number  of 
possible  basic  solutions  can  be  very  large.  It  is  therefore  desirable  to  quickly  find  the  basic  feasible  solution 
which  minimises  the  objective  function  without  examining  all  possibilities. 

The  basic  idea  of  the  simplex  method  is  to  start  with  a  basic  feasible  solution  and  then  seek  a  neighbouring 
basic  feasible  solution  which  reduce  the  value  of  the  objective  function  further  [40]. 

The  optimisation  procedure  is  a  computational  loop  exploiting  the  (consistently  derived)  analytical  gradient 
and  sensitivities  of  the  objective  function  to  the  control  parameters  (activation  deformation  of  the 
reinforcement  bar  elements)  to  execute  a  series  of  non-linear  finite  element  analyses.  When  the  target  shape 
is  reached  the  optimisation  procedure  is  stopped. 

Remark  3 


The  overall  robustness  of  the  global  procedure  was  found  to  be  mainly  conditioned  by  the  non-linear 
finite  element  solving  procedure.  For  this  reason  a  special  arc-length  solver  capable  of  handling  critical 
points  of  the  equilibrium  path  such  as  limit  points,  bifurcations  points,  snap-back  and  snap-though 
[36],  This  was  detailed  in  section  2.3. 
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2.5  WP4:  Development  and  implementation  of  advanced  algorithms  for  tensor  visualisation 

In  this  work  package,  advanced  eigenvalue/eigenvector  extraction  techniques  and  tensor  visualisation 
algorithms  were  extended/implemented  in  ParaView.  Tensor  field  visualisation  in  solid  mechanics  has  not 
been  fully  embraced  by  researchers  and  engineers  despite  the  physical  insight  that  such  technique  could 
provide  to  interpret  complex  multi-dimensional  data  sets  [41,  42].  Typical  visualisation  techniques  for  tensor 
field  data  consist  of  scalar  surface  colour  plot  of  one  component  of  the  tensor  or  one  of  its  eigenvalues.  In  that 
case,  the  directional  information  contained  in  the  tensor  is  lost.  Typical  tensors  encountered  in  continuum 
mechanics  are  strain,  stress,  permeability,  diffusion,  fabric  or  damage  tensors.  Anisotropic  behaviour,  whether 
it  concerns  mechanical  properties,  residual  stress,  electrical/magnetic  permittivity  can  be  captured  and 
visualised  using  tensor  formalism. 

There  are  different  methods  to  visualise  tensors.  For  example,  one  can  calculate  the  eigenvalues  and 
associated  eigenvectors  of  the  tensor.  Each  eigenvector  defined  at  each  point  of  space  (each  integration  point 
or  each  node  in  the  context  of  finite  element  analyses)  represents  a  vector  field.  Similar  to  classical  flow 
visualisation  techniques,  one  can  integrate  and  plot  lines  tangent  to  the  vector  field.  These  lines  can  be  further 
transformed  into  tubes  (or  other  axially-symmetric  primitives)  to  enhance  the  perception  of  depth  (Figure  4). 
Tubes  can  be  scaled  by  the  associated  eigenvalues.  Another  way  to  represent  a  tensor  field  consists  in  plotting 
a  2D  or  3D  icon  at  each  point  which  accounts  for  the  local  characteristics  of  the  tensor  field.  For  example, 
ellipsoids  are  widely  used  to  represent  symmetric  tensors.  The  principal  axes  of  the  ellipsoid  are  aligned  with 
the  three  eigenvector  directions  and  scaled  according  to  the  corresponding  eigenvalues  (Figure  5). 


Figure  4.  Visualisation  of  the  first 
principal  Green-Lagrange  strain 
vector  field  (using  streamlines) 
and  the  associated  eigenvalue 
(first  principal  strain)  in  a  wrinkled 
isotropic  thin  structure  after 
application  of  an  in-plane 
compression. 


Figure  5.  Close-up  view  of  tensor 
ellipsoids  and  tubed  streamlines 
representing  the  three  principal 
Green-Lagrange  strain  vector  fields 
in  the  structure  presented  in  Figure 
4. 


For  this  project,  the  ability  to  visualise  strain  and  stress  tensors  will  prove  very  useful  when  considering 
complex  arrangements  of  active  elements  and/or  anisotropic  materials  and/or  structures  made  of  pre-stressed 
polymer  sheets  and/or  multiphysics  materials  (electro-  and/or  magneto-active).  Moreover,  novel  techniques 
to  visualise  non-symmetric  tensors  [43]  could  be  also  implemented  in  the  future. 
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3.  Exploitation  of  the  computational  platform 


The  exploitation  of  the  computational  platform  developed  in  the  first  part  of  this  project  (see  section  2) 

demonstrated  the  feasibility  of  the  proposed  concept  for  tunable  shape-shifting  or  morphing  structures. 


3.1  Computational  platform  to  design  a  shape-shifting  or  morphing  structure 

A  typical  computational  workflow  for  determining  surface  tuning  parameters  can  be  broken  down  in  the 
following  steps: 


1.  Build  a  multi-layer  cuboid  solid 


Figure  6.  Solid  representation  of  a 
3-layer  structure  of  dimensions 
axbx(hl+h2+t)  reinforced  by  4  bar 
elements  (BARI,  BAR2,  BAR3  and 
BAR4) 

Bl:  solid  layer  1 
B2:  solid  layer  2 
SI:  shell  layer  1 


2.  Embed  reinforcement  bar  elements  (number  of  bar  reinforcements  can  be  arbitrarily  chosen)  can  be 
in  the  shell  surface  (SI)  according  to  a  specific  geometric  pattern. 

3.  Assign  material  properties  to  each  of  the  layers.  In  this  project,  simple  isotropic  hyperelastic  materials 
(neo-Hookean  formulation  [44])  were  used  but  other  types  of  materials  (e.g.  anisotropic 
hyperelastic/hyperviscoplastic)  could  be  seamlessly  used  in  this  framework.  For  example,  in  any 
follow-up  study  of  this  project,  constitutive  models  for  electroactive  polymers  could  be  implemented. 

4.  Apply  periodic  boundary  conditions  to  the  lateral  faces  of  the  cuboid  and  prevent  vertical 
displacement  of  the  bottom  face  of  the  base  layer  (Bl). 

5.  Specify  a  target  shape  for  the  top  surface  of  the  multi-layer  structure  by  providing  an  analytical 
equation  z  =  f(x,y ) . 

6.  If  needed,  specify  additional  constraints  (e.g.  maximum  extension/contraction  of  reinforcement  bar 
elements). 

7.  Run  the  optimisation  procedure  to  determine  activation  patterns  for  the  reinforcement  bar  elements 
that  lead  to  deformations  matching  the  target  3D  topography  of  the  surface.  The  optimisation 
procedure  is  a  loop  exploiting  the  gradient  and  sensitivities  of  the  objective  function  to  the  control 
parameters  (activation  deformation  of  the  reinforcement  bar  elements)  to  execute  a  series  of  non¬ 
linear  finite  element  analyses.  When  the  target  shape  is  reached  the  optimisation  procedure  is 
stopped 

8.  Obtain  the  activation  patterns  for  the  reinforcement  bar  elements. 


Remark  4 


Computational  cost:  for  the  examples  considered,  the  total  computational  time  of  the  optimisation 
procedure  varied  between  20  and  65  minutes.  No  doubt,  this  could  be  accelerated  using  a  larger  number 
of  CPUS  and/or  exploiting  the  excellent  capabilities  of  modern  GPU  architecture. 
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3.2  Examples  of  morphed  surfaces 

In  this  section  examples  of  morphed  surfaces  are  highlighted.  The  reference  multi-layer  structure  is  the  one 
depicted  on  Figure  6,  reproduced  for  convenience  on  Figure  7.  It  is  important  to  note  that  the  following 
characteristics  were  considered: 

■  The  3  layers  and  active  elements  were  made  of  isotropic  hyperelastic  materials  (neo-Hookean 
formulation) 

■  Solid  layer  Bl:  10  MPa  initial  Young's  modulus(0.45  Poisson's  ratio)  /  1.6  mm  thickness 

■  Solid  layer  B2:  50  MPa  initial  Young's  modulus(0.45  Poisson's  ratio)  /  0.4  mm  thickness 

■  Shell  layer  SI:  500  MPa  initial  Young's  modulus(0.45  Poisson's  ratio)  /  0.08  mm  thickness 


Figure  7.  a/  Multi-layer 
structure  reinforced  by  4 
active  bar  elements  in  its 
reference  configuration; 
b/  Same  structure  after 
active  deformation  of  the 
bar  elements.  The  colour 
corresponds  to  the  vertical 
displacement. 


Active  element:  5000  MPa  initial  Young's  modulus(0.3  Poisson's  ratio) 


b/ 


Figure  8.  Examples  of  target  surface  topographies  obtained  by  optimisation  of  active  deformations  of  reinforcement  bar 
elements.  The  initial  configuration  is  that  showed  in  Figure  7-a.  The  grey  surface  represents  the  top  surface  of  the  multi¬ 
layer  structure  before  application  of  active  deformations. 
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4.  Conclusions  and  opportunities  for  further  research  and  development 


The  modelling  project  supported  by  EOARD  [FA8655-12-1-2103]  demonstrated  that  it  is  possible  to  three- 
dimensionally  deform  the  surface  of  a  multi-layer  soft  structure  to  match  a  given  target  surface  topography  by 
inducing  controlled  localised  buckling.  This  controlled  buckling  is  achieved  by  active  contraction  and/or 
extension  of  active  reinforcement  bar  elements  embedded  in  the  top  surface  layer  of  a  given  polymeric 
material.  Combining  several  layers  (with  distinct  properties)  in  the  same  structure  tremendously  increase  the 
number  of  potential  shape  configurations.  The  computational  platform  that  was  developed  is  general  and  can 
accommodate  various  types  of  active  materials  acting  as  controllable  actuators.  It  is  only  a  question  of 
implementing  the  appropriate  multiphysics  constitutive  models  (e.g.  electro-mechanics)  to  virtually  design 
active  surfaces  made  of  commercially  available  or  purpose-built  electroactive  polymers. 

4.1  Some  thoughts 

At  this  preliminary  stage,  we  have  only  scratched  the  surface  of  what  is  possible.  From  section  3.2  (Figure  8), 

one  can  appreciate  the  variety  of  morphed  surfaces  that  can  be  obtained  with  a  very  simple  arrangement  of  4 
reinforcement  bar  elements. 

It  is  noteworthy  that  what  differentiates  the  obtained  morphed  surfaces  is  solely  the  pattern  of  active 
deformations  in  the  reinforcement  bar  elements  that  was  determined  through  an  integrated  optimisation 
procedure. 

It  is  therefore  legitimate  to  speculate  that  variations  in  individual  layer  thickness  and  mechanical  properties 
(e.g.  increased/reduced  stiffness  or  introduction  of  anisotropic  materials)  could  lead  to  an  even  greater  level  of 
control  for  the  type  of  active  surfaces  considered  in  this  project. 

Also,  during  the  iterative  optimisation  procedure,  before  numerical  convergence  to  the  target  shape, 
intermediate  3D  shapes  are  obtained.  These  shapes  could  be  saved  and  automatically  analysed  (using  for 
example  Fourier  transform  or  multiphysics-finite  element  simulations)  to  extract  potential  interesting  physical 
properties  (e.g.  acoustic,  optical,  mechanical)  that  were  not  anticipated. 

The  computational  platform  can  be  used  to  study  morphing  structures  operating  from  the  nanoscopic  to  the 
macroscopic  scale.  Naturally,  if  considering  small  scales,  surface  physics  phenomena  such  as  van  der  Waals 
forces  [45]  would  have  to  be  included.  Because  of  the  modular  nature  of  the  computational  platform  it  would 
be  straightforward  to  implement  these  types  of  force  using  dedicated  field  variables  besides  displacement. 

4.2  Applications 

As  already  highlighted,  reconfigurable  active  surfaces  have  a  huge  potential  in  a  military  (and  also  civil) 
context.  The  proposed  innovation  could  not  only  accelerate  research  and  development  in  metamaterials  and 
active  surfaces  but  could  also  improve  the  performance  of  existing/future  systems  by  integrating 
reconfigurable  multi-functionality.  Selected  potential  applications  ( all  with  the  added  benefit  of  being 
reconfigurable)  of  the  proposed  research  are: 

■  smart  surfaces  that  can  change  their  radar,  acoustic,  thermal  and/or  optical  signatures  (camouflage) 

■  adjustable  inter-locking  mechanisms  for  under-sea/space  structures; 

■  directional  textured  surfaces; 

■  directional  friction  surfaces  to  enhance  manoeuvrability  of  aerospace  and  under-sea  vehicles; 

■  patterned  surfaces  for  drag  reduction; 

■  anti-biofouling  patterned  surfaces  for  over-  and  under-sea  vehicles/structures 

■  active  cooling 

■  context-sensitive  reconfigurable  control  surfaces  (e.g.  keyboards  switching  from  one  language  to 
another). 
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4.3  Future  work 

This  project  mainly  focused  on  developing  a  numerical  toolbox  and  environment  to  model  reconfigurable 
active  structures.  Exploitation  of  the  modelling  platform  has  been  initiated  but  a  vast  amount  of  exploratory 
work  remains  to  be  done.  Large  scale  parametric  analyses  combining  active  deformations  with  variations  in 
geometrical  and  mechanical  properties  and  also  arrangement  of  active  reinforcement  elements  (e.g.  circular, 
concentric,  butterfly  wings,  spiderweb)  should  be  conducted.  Multiobjective  design  optimisation  could  also  be 
conducted  to  design  surface  topographies  fulfilling  particular  physical  properties  (e.g.  wave  refraction). 

In  this  project  the  physics  triggering  active  deformations  within  the  reinforcement  bar  elements  was  not 
accounted  for.  It  would  therefore  make  sense  that,  in  a  sequel  to  this  project,  realistic  physical  mechanisms  for 
active  deformations  be  implemented  in  the  computational  platform.  One  could  envisage  the  experimental 
characterisation  of  various  formulations  of  electroactive  polymers  to  guide  the  development  of  a  physically- 
sound  multiphysics  constitutive  model  of  electroelasticity.  The  active  deformations  of  simple  samples  made  of 
these  polymers  could  be  tested  and  used  to  validate  the  constitutive  modelling  aspects.  After  validation  of  the 
electroactive  characteristics  of  the  polymers,  one  could  develop  a  physical  multi-layer  prototype  structure  and 
use  the  computational  platform  to  establish  active  control  patterns  necessary  to  obtain  target  shapes.  The 
results  of  the  numerical  optimisation  procedure  could  then  be  used  to  control  the  physical  prototype  and, 
hopefully,  validate  the  design  platform. 
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6.  Appendix  A  -  Journal  paper  submitted 


Paper  submitted  to  the  journal  Computer  Methods  in  Applied  Mechanics  and  Engineering: 

Chen,  L.,  Nguyen-Thanh,  N.,  Nguyen-Xuan,  H.,  Rabzcuk,  T.,  Bordas  S.  P.  A.  and  Limbert,  G.  (2013)  Explicit  finite 
deformation  of  isogeometric  membranes,  Computer  Methods  in  Applied  Mechanics  and  Engineering  (under 
review) 


14 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


Elsevier  Editorial  System(tm)  for  Computer  Methods  in  Applied  Mechanics  and 


Engineering 

Manuscript  Draft 


Manuscript  Number:  CMAME-D-13-00626R1 

Title:  Explicit  finite  deformation  analysis  of  isogeometric  membranes 
Article  Type:  Research  Paper 

Keywords:  Membrane;  Kirchhoff-Love  shell;  isogeometric;  NURBS;  explicit;  dynamic  relaxation 
Corresponding  Author:  Dr.  Georges  Limbert,  MSc,  PhD 
Corresponding  Author's  Institution:  University  of  Southampton 
First  Author:  Lei  Chen 

Order  of  Authors:  Lei  Chen;  Nhon  Nguyen-Thanh;  Hung  Nguyen-Xuan;  Timon  Rabczuk;  Stephane 
Pierre  Alain  Bordas;  Georges  Limbert,  MSc,  PhD 

Abstract:  NURBS-based  isogeometric  analysis  was  first  extended  to  thin  shell/membrane  structures 
which  allows  for  finite  membrane  stretching  as  well  as  large  deflection  and  bending  strain.  The 
assumed  non-linear  kinematics  employs  the  Kirchhoff-Love  shell  theory  to  describe  the  mechanical 
behaviour  of  thin  to  ultra-thin  structures.  The  displacement  fields  are  interpolated  from  the 
displacements  of  control  points  only,  and  no  rotational  degrees  of  freedom  are  used  at  control  points. 
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